#include <string>
#include <algorithm>
#include <iostream>
#include <vector>
#include <sstream>
#include <stdexcept>


// boost functions
#include <boost/lexical_cast.hpp>

// my functions
#include <portable_timer.h>
#include <progress_indicator.h>
#include <postgres/pg_copy_guard.h>
#include <postgres/pg_db_cursor.h>
#include <bl2seq_mine.h>
#include <print_error.h>

#include "get_arg.h"

#include <gc.h>
using std::cerr;
using std::cout;
using std::string;
using std::vector;
using std::sort;
using boost::lexical_cast;
using std::ostringstream;


struct t_gi_seq
{
	unsigned gene_position;
	string gene_id;
	string acc;
	string pi;
	string seq;
	t_gi_seq(unsigned pos, const string& g, const string& a, const string& p, const string& s):
				gene_position(pos), gene_id(g), acc(a), pi(p), seq(s)
	{
		bl_prepare_sequence_str(seq);
	}
};

typedef vector<vector<t_gi_seq> > t_vec_vec_gi_seq;

void save_alignment_to_db(	t_pg_copy_guard&		to_db,
							const t_gi_seq&			prot1,
							const t_gi_seq&			prot2,
							t_blast_results&		blast_res);

//________________________________________________________________________________________

//		align

//________________________________________________________________________________________
void align(t_pg_copy_guard& to_db, const t_vec_vec_gi_seq& vec_vec_gi_seq,
					unsigned gene1, unsigned gene2,
					t_blast_options& blast_options_seg, t_blast_options& blast_options,
					unsigned& cnt_blast_comparisons)
{
	unsigned best_i = 99999;
	unsigned best_j = 99999;

	//
	//	Find best 'seg'ed alignment
	//
	t_blast_results best_blast_res;
	for (unsigned i = 0; i < vec_vec_gi_seq[gene1].size(); ++i)
	{
		for (unsigned j = 0; j < vec_vec_gi_seq[gene2].size(); ++j)
		{
			++cnt_blast_comparisons;
			//
			//	blast;
			//
			t_blast_results blast_res;


				if (!do_bl2seq(	blast_options_seg,
							vec_vec_gi_seq[gene1][i].seq,
							vec_vec_gi_seq[gene2][j].seq,
							blast_res)				||
					blast_res.align_len < 50)
				continue;

			//
			//	Is this a better alignment?
			//

			// save only if bitscore larger
			if (blast_res.bitscore < best_blast_res.bitscore)
				continue;
			// if same evalue
			else if (blast_res.bitscore == best_blast_res.bitscore)
			{
				//only save if lower e_value
				if (blast_res.e__value > best_blast_res.e__value)
					continue;

				//  or shorter length
				if (blast_res.total_len  >= best_blast_res.total_len)
					continue;
			}

			//
			//	Save better alignment
			//
			best_blast_res 	= blast_res;
			best_i			= i;
			best_j			= j;
		}
	}

	if (best_i == 99999)
		return;



	//
	//	Find best un'seg'ed alignment
	//
	t_blast_results blast_res;
	if (!do_bl2seq(	blast_options,
					vec_vec_gi_seq[gene1][best_i].seq,
					vec_vec_gi_seq[gene2][best_j].seq,
					blast_res))
		return;
	save_alignment_to_db(	to_db,
							vec_vec_gi_seq[gene1][best_i],
							vec_vec_gi_seq[gene2][best_j],
							blast_res);

}


//________________________________________________________________________________________

//		save_alignment_to_db

//________________________________________________________________________________________
void save_alignment_to_db(	t_pg_copy_guard&		to_db,
							const t_gi_seq&			prot1,
							const t_gi_seq&			prot2,
							t_blast_results&		blast_res)
{
//	"(align_start, align_len, aligned_positions, coverage, "
//	"e_value, len1, len2, aligned_residues, bitscore, "
//	"gaps1, gaps2)\n\n";
	ostringstream ostrm;
	if (prot1.gene_id < prot2.gene_id)
		ostrm
			<< prot1.gene_id                        << "\t"
			<< prot2.gene_id                        << "\t"
			<< prot1.acc                            << "\t"
			<< prot2.acc                            << "\t"
			<< prot1.pi                             << "\t"
			<< prot2.pi                             << "\t"
			<< prot1.gene_position                  << "\t"
			<< abs(	prot1.gene_position -
					prot2.gene_position)            << "\t"
			<< blast_res							<< "\n";
	else
	{
		// the other way around
		blast_res.swap_subject_query();
		ostrm
			<< prot2.gene_id                        << "\t"
			<< prot1.gene_id                        << "\t"
			<< prot2.acc                            << "\t"
			<< prot1.acc                            << "\t"
			<< prot2.pi                             << "\t"
			<< prot1.pi                             << "\t"
			<< prot2.gene_position                  << "\t"
			<< abs(	prot2.gene_position -
					prot1.gene_position)            << "\t"
			<< blast_res							<< "\n";
		// swap back just in case
		blast_res.swap_subject_query();
	}

	to_db += ostrm.str();
	to_db.cached_write();
}











//________________________________________________________________________________________

//		get_gene_sequences

//________________________________________________________________________________________
bool get_gene_sequences(t_pg_cursor& pg_db, const string& species,
						const string& chromosome,
						t_vec_vec_gi_seq& vec_vec_gi_seq,
						t_program_args& args)
{
	string sql_cmd = "SELECT gene_position, gene_id, acc_code, pi, sequence "
						"FROM "
							"sequences NATURAL JOIN "
							"ens_id NATURAL JOIN "
							+ species  + ".gene_loci "
						"WHERE "
							"chromosome = '" + chromosome + "' "
						"ORDER BY gene_position";
	// create binary cursor
	if (!pg_db.Declare(sql_cmd, false))
		return pg_print_sql_err("create a cursor to "
								"retrieve sequences for chromosome " + chromosome,
								pg_db,
								sql_cmd, args.CERR());

	{
		unsigned gene_position = 5000000;
		unsigned batch_size = 500;
		// get data in batches
		while (pg_db.Fetch(batch_size))
		{
			int curr_batch_size = pg_db.Tuples();
			if (!curr_batch_size)
				break;
			for (int i = 0; i < curr_batch_size; ++i)
			{
				unsigned new_gene_position = lexical_cast<unsigned>(pg_db.getvalue(i, 0));
				if (new_gene_position !=gene_position)
				{
					vec_vec_gi_seq.push_back(vector<t_gi_seq>());
					gene_position = new_gene_position;
				}
				vec_vec_gi_seq.back().push_back(t_gi_seq(   gene_position,
															pg_db.getvalue(i, 1),
															pg_db.getvalue(i, 2),
															pg_db.getvalue(i, 3),
															pg_db.getvalue(i, 4)));
			}

		}

	}

	// close cursor
	if (!pg_db.Close())
		return pg_print_sql_err("close the  cursor for "
								"retrieve sequences for chromosome " + chromosome,
								pg_db, "CLOSE CURSOR", args.CERR());
	return true;
}

//________________________________________________________________________________________

//		get_chromosomes

//________________________________________________________________________________________
struct cmp_chromosomes
{
	bool operator ()(string a,  string b) const
	{
		if (a.length() < b.length())
			a = string(b.length() - a.length(), ' ') + a;
		else if (b.length() < a.length())
			b = string(a.length() - b.length(), ' ') + b;
		return a < b;
	}
};

bool get_chromosomes(t_pg_cursor& pg_db, const string& species,
						vector<string>& chromosomes,
						t_program_args& args)
{
	bool all_chromosomes = args.chromosome == "all";
	string sql_cmd = "SELECT distinct on (chromosome) chromosome "
						"FROM "
							+ species  + ".gene_loci "
						"WHERE chromosome "	+
						(all_chromosomes 				?
							string("~ '^[0-9XY]+$'")	:
   							("= '" + args.chromosome + "'")
						 ) +
						" ORDER BY chromosome";
	if (!pg_db.ExecTuplesOk(sql_cmd))
	{
		pg_print_sql_err("retrieve the list of chromosomes for " + species +
							" from postgreSQL",
							pg_db, sql_cmd, args.CERR());
		return false;
	}

	unsigned cnt = pg_db.Tuples();
	if (!cnt)
	{
		throw std::logic_error("No matching chromosomes for " + species + " were found in the panda "
				+ (all_chromosomes ? string() : ( " which matched '" + args.chromosome + "'")));
	}
	for (unsigned i = 0; i < cnt; ++i)
		chromosomes.push_back(pg_db.GetValue(i, 0));
	sort(chromosomes.begin(), chromosomes.end(), cmp_chromosomes());
	return true;
}

void do_nothing_for_gc_warnings(char *msg, GC_word arg)
{
}


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		Main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main(int argc, char* argv[])
{
	GC_set_warn_proc(do_nothing_for_gc_warnings);
	t_portable_timer timer;

	// don't use C streams
	std::ios::sync_with_stdio(false);

	// saving parameters
	string exe_parameters;
	for (int i = 0; i < argc; ++i)
		exe_parameters += string(argv[i]) + ' ';

	// we use cout to print progress so flush per call
	cout.setf(std::ios_base::unitbuf);

	// get command line arguments
	t_program_args args;
	try
	{
	if (!process_command_line_options(argc, argv, args))
		return 1;


	//____________________________________________________________________________________
	// 		Connect to PostgreSQL
	//____________________________________________________________________________________
	t_pg_cursor pg_db(args.db_login.c_str(), "panda_cursor");	//"dbname=all_vs_all_gene_bl2seq user=postgres"
	if ( pg_db.ConnectionBad() )
	{
		cerr	<< "\nERROR!! Connection to database ["
							<< args.db_login << "] failed.\n"
							<< pg_db.ErrorMessage() << "\n\n";
		return 1;
	}

	//____________________________________________________________________________________
	// 		get largest cluster id
	//____________________________________________________________________________________


	vector<string> chromosomes;
	if (!get_chromosomes(pg_db, args.species, chromosomes, args))
		return 1;

	t_blast_options blast_options_seg(true, 1e-5, "BLOSUM80");
	t_blast_options blast_options(false, 1e-5, "BLOSUM80");
	unsigned cnt_blast_comparisons = 0;
	unsigned cnt_genes = 0;
	// write per chromosome
	for (unsigned i = 0; i < chromosomes.size(); ++i)
	{
		t_vec_vec_gi_seq vec_vec_gi_seq;
		args.VERBOSE() << "Chromosome " + chromosomes[i] + ": ";
		if (!get_gene_sequences(pg_db, args.species, chromosomes[i], vec_vec_gi_seq, args))
			return 1;
		cnt_genes += vec_vec_gi_seq.size();


		t_pg_copy_guard copy_buffer(pg_db, args.output_table +
			"(gene_id1,gene_id2, acc_code1, acc_code2, pi1, pi2,"
									" gene_position, gene_loci_diff, " +
									t_blast_results::get_sql_column_names() + ")",
									args.CERR());
		if (!copy_buffer)
			return 1;

		args.VERBOSE() << "Aligning " << vec_vec_gi_seq.size() << " adjacent genes ...\n";
		t_progress_indicator dots(args.VERBOSE(), 5, vec_vec_gi_seq.size());
		dots.use_timer();
		for (unsigned curr_pos = 0; curr_pos + args.window_size < vec_vec_gi_seq.size(); ++curr_pos, ++dots)
			for (unsigned offset = 1; offset < args.window_size; ++offset)
				align(copy_buffer, vec_vec_gi_seq, curr_pos, curr_pos + offset,
						blast_options_seg, blast_options, cnt_blast_comparisons);

		if (!copy_buffer.end_copy())
			return 1;
		pg_db.Commit();
		pg_db.BeginTransaction();
	}


	args.VERBOSE() << string(80, '_') << "\n";
	args.VERBOSE() << cnt_genes << " genes in total\n";
	args.VERBOSE() << timer.reset() << " Seconds for " << cnt_blast_comparisons << " blast comparisons\n";
	}
	catch (std::exception& e)
	{
		std_print_error(args.CERR(), e.what());
		return 1;
	}

	return 0;
}



